namespace IBS.Math

module Solver =
    open System

    let private precision = 1e-7

    // a is min bound, b - max
    let private findBounds f guess a b =
        if guess <= a || guess >= b then failwith (sprintf "Guess should be inside the interval (%f, %f)" a b) 
        let shift, k = 0.01, 1.6

        let inline (!<) v = if v <= a then a + precision else v
        let inline (!>) v = if v >= b then b - precision else v
        
        let rec find low up = function
            | 0     -> failwith (sprintf "The bounds can't be found - the last interval (%f, %f)" low up)
            | iters ->
                 let low, up = !< low, !> up
                 match f low, f up with
                 | 0., _ | _, 0.      -> low, up
                 | x, y when x*y < 0. -> low, up
                 | x, y when x*y > 0. -> find (low + k * (low - up)) (up + k * (up - low)) (iters - 1) 
                 | _                  -> failwith "Function result for one of the limits isn't a plain number"
        
        let low = !< (guess - shift) 
        let high = !> (guess + shift)        
        find low high 50

    let private d f x = (f (x + precision) - f (x - precision))/(2. * precision)
    
    let rec newton f x iters =  
        let fx, Fx = f x, d f x
        let x' = x - (fx / Fx)
        if abs (x' - x) < precision then Some x'
        elif iters > 50 then None
        else newton f x' (iters + 1)

    let rec bisection f a b count =
        if a = b then failwith "Invalid bisection interval" 
            
        let fa = f a
        if abs fa < precision then a // a is the root
        else
            let fb = f b
            if abs fb < precision then b // b is the root
            else
                let newCount = count + 1
                if newCount > 20 || fa * fb > 0. then failwith "Bisection: the root can't be found"

                let mid = a + 0.5 * (b - a)
                let fmid = f mid
                    
                if abs fmid < precision then mid // the root
                elif fa * fmid < 0. then bisection f a mid (count + 1)
                elif fa * fmid > 0. then bisection f mid b (count + 1)
                else failwith "Bisection: the root can't be found" 

    let zeroin (ax, bx) f tol =
        let eps = 1e-5
        // swap data for b to be the best approximation
        let swap a b c =
            let fa, fb, fc = f a, f b, f c
            if abs fc < abs fb then
                b, c, b, fb, fc, fb
            else 
                a, b, c, fa, fb, fc
                
        let inline checkSigns (q, p) = if p > 0. then -q, p else q, -p        
        let inline adjustStep s t = if abs s < t then float (sign s) * t else s
                    
        let rec iter a b c =
            let prevstep = b - a //distance from the last but one to the last approximation        
            let a, b, c, fa, fb, fc = swap a b c
            let actTolerance = 2. * eps * abs b + tol / 2. // actual tolerance
            let c_b = c - b
            let newstep = c_b / 2.
            
            if abs newstep <= actTolerance || fb = 0. then
                b //acceptable approximation
            else
                //decide if interpolation can be tried - prevstep was large enough and was in true direction
                let step = 
                    if abs prevstep < actTolerance || abs fa <= abs fb then newstep
                    else
                        let fba = fb / fa
                        let q, p =
                            if a = c then
                                1. - fba, c_b * fba
                            else
                                let fac, fbc = fa / fc, fb / fc
                                (fac - 1.)*(fbc - 1.)*(fba - 1.), fba*(c_b*fac*(fac - fbc) - (b - a)*(fbc - 1.))
                            |> checkSigns

                        // if b + p / q falls in [b, c] and isn't too large it's accepted
                        // if p / q is too large then the bissection procedure can reduce [b,c] range to more extent
                        if p < 0.75 * c_b * q - abs (actTolerance * q) / 2. && p <  abs (prevstep * q / 2.) then
                            p / q
                        else newstep
                
                // adjust the step to be not less
                // save the prev approximation
                let b' = b + adjustStep step actTolerance
                let c = if sign (f b') = sign fc then b else c
                iter b b' c
        iter ax bx ax   
    
    // use combination of bisection and interpolation to find a root
    let findRoot f guess =
        match newton f guess 0 with
        | Some r -> r
        | _ -> zeroin (0., 3.) f 1e-7 

